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ABSTRACT | 
fa Naa ICWAKE computer code solves the Navier-Stokes equations for 


| axisymmetric, incompressible, swirling, turbulent flow with large «xial 
gradients. This document is a guide to the use of the code. Included 
are descriptions of the input parameters and the code structure, sone 


8 general comments about using the code and a sample calculation. 
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1. INTRODUCTION 

This manual describes the structure of the ICWAKE computer code 
and the use of this code to solve the Navier-Stokes equations for 
axisymmetric, incompressible, swirling, turbulent flow with large 
axial gradients. The mathematical analysis and finite-difference 
formulation upon which it is based are described in the companion 
document "ICWAKE COMPUTER CODE —— Mathematical Analysis and Finite- 
Difference Formulation" (Grabowski, et al. 1976). Prospective users of 
the code should be familiar with this reference. Experience in the 
use of large hydrodynamics computer codes is necessarily assumed be- 
cause the code obtains a finite-difference solution for an elliptic 
system of eleven rather lengthy, coupled, partial differential equa- 
tions. 

ICWAKE, which is written in FORTRAN IV, was developed and used 
extensively on the Lawrence Berkeley Laboratory CDC 7600 computing 
system. It can be compiled to produce executable object codes on 
the LBL RUN76, FIN and FIN4 compilers. ICWAKE requires approximately 
160K storage to execute with a maximum grid system of 60 x 32 points. 
It makes use of the NCAR direct Poisson solver, BLKTRI. 

ICWAKE is based on a code, FIND@OM, for laminar flow, which was 
written under the sponsorship of NASA/Ames Grant NGR 05-003-45! at 
the University of California, Berkeley. Further development ard the 
extension to turbulent flow were performed under the sponsorship of 
a Flow Research, Inc., Independent Research and Development Program, 
and most recently under the sponsorship of ONR Contract No. NO0014- 
76-C-0564. 

In section 2 of this manual, we describe the required input 
quantities for code execution. In section 3, we outline the struc- 
ture of the code and the roles of the main routine and the subroutines. 
In section 4, we present some general observations about the use of 
the code and the choice of parameters. Finally, in section 5 we present 


input and cutput for a complete calculation. 
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2. DESCRIPTION OF INPUT QUANTITIES 
‘ The input data necessary to exercise ICWAKE enter the code through 
Namelists DAT1 through DAT6, DIS7, and DAT8, which are read by the main 
program, subroutine DATIN, and subroutine STRESS2, respectively. All 
input parameters for which no dimensions are specified are dimensionless. 
. Data from previous calculations can be entered through logical unit 41 
(TAPE41), which is read by ICWAKE, and can be used as the starting con- 
ditions. (Using data from previous calculations as initial conditions 
is discussed in the description of routines UPC@ND, DATIN, and IUNI in 
c Section 3.) 
Namelist $DAT1 
$DAT1 is read at the beginning of each exercise of the code. It 
provides the following: the parameters that determine the frequency and 
: scope of the printout; the time-step size, and a maximum number of time- 
steps to be taken in the computation; a convergence criterion, and the 1 
frequency at which convergence checks are performed; the sequence number 
of the computer run, and the time-step count at the beginning of the 
. run. 
$DAT1 
M1,M2,M3 Output printout control at intermediate time-steps. { 
Solution is printed out for all I (i.e., x) at values of 
c J from Ml through M2 in increments of M3 (Default values 
= 1,1,1) i 
M11,M22,M33 Output printout control at termination of the run (Default 
values = 1,8,1) 
c NPRNT Frequency (in terms of time-steps) of intermediate 
printout (Default value = 30) 
TA Time-step size (Default value = 0.1) 
NTMX Maximum number of time-steps to be taken (Default value = 30) 1 
CONCRIT Convergence criterion based on the time rate of change of 
the rms divergence of mean velocity (Default value = 1.E-9) 
NTCHK Frequency (in terms of time-steps) of convergence checks 
(Default value = 10) 
NUMBER Sequence number of the run (Default value = 1) 
NSTRT Time-step count at the beginning of the run (Default | 
value = 0) 
‘ 
a 7 | 
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Namelist $DAT2 
$DAT2 is read at the beginning of each exercise of the code in 


which data from previous calculations is not used as input (i.e., 


TAPE41 not read). It contains the parameters necessary to set up 


the computational grid system and to generate the coefficients of the 


logarithmic transformation. Generating these four coefficients re- 


quires specifying the total number of grid points desired in each 


direction, as well as the number of points desired between ZINITL 
and some distance ZCC (less than ZMAX), and between the axis of 


symmetry and some distance RCC (less than RMAX). 


t 
{ 
{ 
t 
SDAT2 
N 
c M 
NC 
ct 
MC 
ZINITL 
, 
ZCC 
F RCC 
ZMAX 
C 
YMAX 
P AX, AY 
EPS 
q 


Total number of grid points in the x (or z) direction 
(Default value = 40) 


Total number of grid points in the y for r) direction. 
M must equal 2", where k is an integer greater than 


one (Default value = 16) 


Number of grid points in the x direction required be-~ 
tween ZINITL and ZCC (Default value = 10) 


Number of grid points in the y direction required be~ 
tween the axis (r = 0) and RCC (Default value = 8) 


The z location of the WEST (upstream) boundary. This 
will correspond to x = 0.(Default value = 0) 


A length in the z direction beginning at z= 0. NC 
points are to be located between ZINITL and ZCC. 
(ZCC must be greater than ZINITL) (Default value = 1) 


A length in the r direction beginning at the axis, in 
which MC points are to be located (Default value = 1) 


z location of the EAST (outflow boundary) (at ZMAX, 
x = 1) (Default value = 20) 


r location of the NORTH (radial) boundary (Default 
value = 6) 


y value assigned to the radial boundary itMAX (Default 
value = 0.5) 


Initial estimates of the axial and radial logarithmic 
transformation parameters (Default values = 3.076 and 
5.8, which correspond to the default grid system) 


Convergence coefficient for the logarithmic transforma- 
tion coefficient calculation procedure (Default value = 1.E-6) 


seep 
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MTURB 
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The NORTH boundary conditions on all turbulence 
quantities are applied at J = MTURB (Default 
value = M) 


Fi ere pnp 


t ame 
Namelist $DAT3 
$DAT3 provides the flow Reynolds number, a laminar-turbulent flag, 
s a flag that may be set to decouple the mean flow calculation from the 
turbulence calculation, and a flag that may be set to apply Mager's (1972) 
cubic-quartic profiles at the WEST boudnary. Note that $DAT3 is not 
read during calculations with NSTRT greater than zero. 
t 
SDAT3 
RE Flow Reynolds number based on the assumed character- 
istic velocity and length 
i IDECGUP A flag that, if set to unity, will cause the mean 
flow calculations to be decoupled from the turbulence 
calculation (Default value = 0) 
ITURB A flag that must be set to either zero or unity for 
laminar or turbulent flow, respectively (No default 
t value) 
IMAGER A flag that is set to unity in order to apply Mager's 
cubic-quartic upstream profiles for a laminar calcula- 
tion (Default value = 0) 
C WI Free-stream axial velocity for Mager's quartic profile 
VI Free-stream circulation for Mager's cubic profile 
ALPH Form factor for Mager's quartic profile 
C Note that WI, VI, and ALPH have to be set only when IMAGER = 1. They have 
no default value~. 
The following parameters are required only when IBLPRP = 1 (see 
namelist $DAT4); they should be given in CGS units. (No default values 
C are assigned.) 
| BRD Body radius 
BLN Body length 
( BUIN Body velocity (free-stream axial velocity) 
RPROP Propeller radius (distance from the axis of revolution 
to the propeller tip) 
ALAMB Advance ratio of the propeller, \ = BUIN/(RPROP*2), 
C where 2 is the angular velocity of the propeller 
i (rad/sec) 
NB Number of blades on the propeller 
¢ 


Flow Research Note No. 95 
November 1978 


a le ee eee ee nee 


NBEL 


BELR 


BELC 


BELANG 


BELTH 


BELDEL 
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Number of blade elements 


An array containing the radius of each blade element 


An array containing the chord length of each blade 
element 


An array containing the geometric pitch angle (in 
degrees) of each blade element 


An array containing the maximum thickness of each blade 
element 


An array containing the maximum displacement of the 
mean camber line from a chord line connecting the lead- 
ing and trailing edges 
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Namelist $DAT4 

$DAT4 provides radial profiles of the mean flow and turbulence 
quantities at z = ZINITL (or correspondingly, x = 0). S$DAT4 is not read 
when IMAGER = 1. If a laminar flow (but not Mager's profiles) is to be 
computed (i.e., IMAGER = 0, ITURB = 0), only MP@INT, RT, UL, VL, and WL 
must be specified. If a turbulent flow is to be computed and IBLPRP = 0 
(see below), the $DAT4 arrays are needed. If some are not available, 


subroutine DATIN should be modified to compute or estimate them. 


SDAT4 
MPQPINT The number of data points at the WEST boundary 
RT WEST boundary data point location array (MP@INT values) 
(Note that RT(1) must equal zero and that when $DAT4 is 
read, RMAX is set to RT(MP@INT). Also, note that ali of 
the arrays below consist of MP@INT elements defined at 
the RT data points). | 
UL Array of radial velocities, U. 
VL Array of circumferential velocities, V. 
WL Array of axial velocities, W. 
UU Array of radial-radial velocity fluctuation correlations, R. 
VV Array of circumferential-circumferential velocity fluc- 
tuation correlations, Rog: 
Ww Array of axial-axial velocity fluctuation correlations, a 
Wy Array of axial-radial velocity fluctuation correlations, Roe 
WV Array of axial-circumferential velocity fluctuation 
correlations, Roe° } 
' 
UV Array of radial-circumferential velocity fluctuation 


| 
F 
a 
fe 
| 


correlations, Rio: 
(No default values are assigned to the above arrays.) 


ISCHETZ When this signal is set at 1, the R.. turbulence correlation 
is initialized according to an eddy viscosity formulation 
(Default value = 0) 


IDEP When this signal is set at 1, the dissipation rate ¢€ is 
initialized by setting it equal to the production rate. 
If IDEP = 0, then ec = Kk? 72, where k = | 
(RL. + Rog + Ro.” and K. and 2 are assigned through 
namelist $DIS7. (Default value = 0) 


emma iam a EI = 
Flow Research Note No. 95 


November 1978 
-8- ' 


IBLPRP When this signal is set at 1, all variables are ini- 
tialized according to the boundary-layer/propeller 
algorithm contained in subroutines BLAYER, PR@PWV, 
PROPU, and PTURB (see namelist $DAT3). (Default 
value = 0) 


ba ] 
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Namelist $DAT5 

$DAT5 provides main program iteration-sweep control. Instead of 
the ADI procedure, either horizontal or vertical line-by-line iteration 
may be specified. (The default ADI procedure is strongly recommended. ) 


$DAT5 

ISWEEPX A flag that should be set to zero to skip horizontal 
ADI sweep (Default value = 1) 

ISWEEPY A flag that should be set to zero to skip vertical 


ADI sweep (Default value = 1) 


Note that the computation will abort if both ISWEEPX and ISWEEPY equal zero. 
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Namelist $DAT6 
$DAT6 provides the axial convection weighting factors as a function 


of axial location. 


$DAT6 


ARTVIS An array (N values) of axial convection weighting 
factors as x (or z) locations. ARTVIS may be set 
to zero, unity, or some function between, for purely 
central, purely upwind, or combined differencing, | 
respectively. ARTVIS should always be set to unity 
at x grid points N and N-1, and should vary smoothly 
.to these values. 
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Namelist $DIS7 


$DIS7 


TSCALE The value 2 in the dissipation rate equation, &€ = Kk?! ? 79 


(see IDEP in namelist $DAT4). (Default value = 0.2) 


XKECON The value Ke in the above equation for the dissipation 
rate (Default value = 0.53) 


ee 
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Namelist $DAT8 

$DAT8 is read by STRESS2. It provides turbulence constants and 
also sets flags that determine the form of the boundary condition 
applied to the circumferential-circumferential velocity fluctuation 
correlation at the axis, in addition to flags that determine the 


particular turbulent diffusion model to be used. 


$DAT8 
NHL A flag set to zero or unity for Daly-Harlow or 
Hanjolic-Launder turbulent diffusion, respectively. 
CEPS Turbulence model constant (Default value = 0.15) 
CEPS1 Turbulence model constant (Default value = 1.44) | 
CEPS2 Turbulence model constant (Default value = 1.90) | 
cso Turbulence model constant (Default value = 0.25) | 
csl Turbulence model constant (Default value = 0.11) 
CSN Turbulence model constant equal to CSO if NHL = 0, | 
and CSN if NHL = 1 
CON1 Turbulence model constant (Default value = 1.5) | 
CON2 Turbulence model constant (Default value = 0.4) 
ISBCTT Boundary condition control at R= 0. Setting 


ISBCTT = 0 requires OR, ,/dr = 0, Rog = R,, at r = 0, 
while setting ISBCTT = I requires 9R,./dr = 
ORog/ or = 0 at r = 0 (Default value = 1) 


a 
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3. ICWAKE CODE STRUCTURE 
3.1 Main Program ICWAKE 

The main program ICWAKE controls the code. It handles nearly 
all input and output, and calls the subroutines that set up the finite- 
difference grid, evaluate the transformation coefficients, calculate 
the effect oi a boundary layer and a propeller, and interpolate, from 
the input data, for the boundary values at the grid points along the 
WEST boundary. ICWAKE reads TAPE41 if data from a previous calculation 
is to be used as the initial conditions. ICWAKE performs the finite- 
difference solution of the mean-flow equations and computes the diver- 
gence of the velocity field; after each half time-step (i.e., after 
each ADI sweep), it calls subroutine STRESS2, which calculates Reynolds 
stresses, and subroutine PRESSR, which computes the pressure. ICWAKE 
calls routines that check for convergence, compute the net mass flux 
into the domain, and compare various terms in the equations of motion. 
At the end of a calculation, ICWAKE writes data on TAPE42, which may 
subsequently be used, through TAPE41, to restart the calculation or 


as initial conditions for a calculation. 
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COMPUTE BOUNDARY LAYER AND PROPELLER 
INFLUENCE AT UPSTREAM BOUNDARY 

BY CALLING BLAYER, PROPWV, PROPU, 

PTURB 


GENERATE GRID SYSTEM, 
COMPUTE TRANSFORMATION COEFFICIENTS, 
CALL COEFF 
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BEGIN ITERATIONS 
DO 666 NT = N1,N2 


T 
ISWEEPY = 0 


VERTICAL SWEEP 
IMPLICIT IN X 
DO 501J = 2,M-1 


COMPUTE TRIDI COEFFICIENTS 
AU ,AV MM, BU, +. -5C0,...DU... 
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COMPUTE REYNOLDS STRESSES 
CALL STRESS2 


.| END OF VERTICAL SWEEP 


HORIZONTAL SWEEP 
IMPLICIT TN Y 
DG 601 I = 2,N 


COMPUTE TRIDI COEFFICIENTS. 
Pe aes.” eae 
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G 
q 
COMPUTE 
DU,DV,DW 
1 
t 
, : 
(s) INVERT TRIDI MATRICES 
r FOR U,V AND W 
t | 
R) 
| 
J 
| 
t | 
------@ 
ia 


es 


COMPUTE REYNOLDS STRESSES 
CALL STRESS2 


Sa aa 


ee 


ss 
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COMPUTE P AT 
EAST BOUNDARY 


COMPUTE P 
CALL PRESSR 


END OF HORIZONTAL SWEEP 


CHECK FOR CONVERGENCE 
CALL DIVC@N 


PRINT SOLUTION IF: 
ICON=1, OR MON(NT,NTRNT)=0, 
OR NT=N2 


COMPUTE NET MASS FLUX 
(IF SOLUTION PRINTED) 
CALL FLXCHK 


ICQN=1 


I 
é 


(TERMINATOR OF TIME-STEP DO-LOOP) 
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SAVE DATA ON TAPE42 


| 
| 
| 
i 


TERM COMPARISONS 
CALL TERMCHK 


: | 
C : on Be: , | ! 
c k 
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SUBROUTINE CQEFF. Subroutine C@EFF calculates the transformation co- | 
efficients that appear in the mean, turbulence model and in the pressure 


equations in transformed form. CQ@EFF calls subroutine NEWT, which 


computes the coefficients of the logarithmic transformations. 


COMPUTE COEFFICIENTS OF 
LOGARITHMIC TRANSFORMATION IN R, 
CALL NEWT 


COMPUTE COEFFICIENTS OF 

LOGARITHMIC TRANSFORMATION IN Z, 
CALL NEWT 

} 


COMPUTE COEFFICIENTS APPEARING 
IN EQUATIONS 


~ 
Sie alate eee eee 


SUBROUTINE NEWT. Subroutine NEWT uses a Newton-Raphson technique to 
compute the coefficients of the logarithmic transformation based on the 
estimated values AX and AY, and on RCC, ZCC, ZINITL, ZMAX, RMAX, N, M, 


> 
2 


os 


NC, and MC. 


eee 
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SUBROUTINES UPCOND, DATIN and IUNI. Subroutine UPC@ND sets the initial 
conditions over the solution domain and computes Q°w/ dre at the WEST 
boundary. If MAGER=1, Mager's cubic-quartic profiles are assumed up- 
stream; otherwise, tabular data read into JWBL is assumed. If a 
solution from a previous calculation is not used as input to a new 
calculation (i.e., NUMBER=1, NSTRT=0), the conditions at the WEST 
boundary are applied everywhere, and subroutine PRESSER is called to 
provide the pressure field. If a calculation with new conditions (i.e., 
new Re or WEST boundary conditions) is begun with data from a previous 
calculation (i.e., NUMBER>1, NSTRT=0) as an initial condition, the 
solution at only the WEST boundary is changed. If a calculation is to 
be continued (i.e., NUMBER>1, NSTRT>O) beginning with existing data 
from TAPE41, only a7 u/lde" is computed. When NSTRT=0 and IMAGER=0 so 
that tabular data is assumed, subroutine DATIN is used to interpolate 
for the WEST boundary values at grid points. 

Subroutine DATIN calls IUNI, which performs either first- or 
second-order interpolation to calculate boundary values at the WEST 
grid points from tabular data. DATIN also estimates upstream conditions 


that are not specified, and prints out the upstream conditions. 
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IMAGER=1 


| | 


COMPUTE UPSTREAM 
CONDITIONS BASED ON 
MAGER'S PROFILES 


COMPUTE 32w/dr2 
AT WEST BOUNDARY 


NUMB 


ER>1 
_NSTRT=0 


NSTRT>O 


BASE INITIAL CONDITIONS 
ON UPSTREAM CONDITIONS 
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COMPUTE DIV 


| 


COMPUTE P 
CALL PRESSR 
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ESTIMATE UNSPECIFIED QUANTITIES 
E.G., EPSL 


‘ INTERPOLATE FOR VALUES AT 
GRID POINTS ALONG WEST BOUNDARY 


COMPUTE, IF NECESSARY, MISSING REYNOLDS 
STRESSES AT WEST BOUNDARY FROM 
EDDY VISCOSITY OR OTHER ASSUMPTION 


PRINT UPSTREAM CONDITIONS 


CSC CORP FIRST- OR SECOND-ORDER 
INTERPOLATION PROGRAM 
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SUBROUTINE STRESS2. Subroutine STRESS2 computes Reynolds stresses 
using the second-order closure turbulence model of Hanjalic and 
Launder (1972) with the pressure/mean-strain correlation of Launder, 


Reece, and Rodi (1975). STRESS2 sweeps in x (implicit in y) only. 


COMPUTE COEFFICIENTS 
ON FIRST PASS 
hota Sleds hee BEGIN SWEEP 
r D¢ 1000 I=2,N 


COMPUTE ELEMENTS OF 
TRIDI MATRICES 


ARR,BRR,CRR,ATT,...,AZZ,.-+5 
Pr a ae. eee 


BEGIN CALCULATION OF 
INHOMOGENEOUS TERMS 
DEP,DRR,DTT,...etc. 


COMPUTE SHEAR STRESS AND DISSIPATION RATE 
PRODUCTION RATES PEP,PRP,PTT,...etc. 


COMPUTE PRESSURE-STRAIN CORRELATION TERMS 
PHIRR,..., and PH2RP... 
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COMPUTE DISSIPATION RATES 


COMPUTE RIGHT-HAND SIDE 
TURBULENT DIFFUSION TERMS 


COMPUTE RIGHT-HAND SIDE 
CONVECTION TERMS 


ASSEMBLE DRR,DTT,...etc. 


APPLY NORTH AND SOUTH 
BOUNDARY CONDITIONS 


INVERT TRIDI MATRICES 
FOR EPSL,TRR,...etc. 


* 
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SUBROUTINE PRESSR. Subroutine PRESSR solves the Poisson equation 


for pressure obtained by taking the divergence of the momentum equa- 
tions. It calls BLKTRI, the NCAR fast direct solver. A flow chart 
for BLKTRI is not presented here. Note that BLKTRI is initialized by 


calling in once with IFLG=0. 


INITIALIZE COEFFICIENTS 
ON FIRST PASS 


COMPUTE LAMINAR INHOMOGENEOUS TERMS 
THESE ARE STORED IN ARRAY P(I,J) 


INCLUDE NORTH, EAST, AND WEST BOUNDARY CONDITIONS 
FOR LAMINAR FLOW IN P(I,J) 
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CALL BLKTRI IFLG=1 


Pd ,J) NOW 
CONTAINS PRESSURE 


Seen eeeneneresere= 
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SUBROUTINE TERMCHK. Subroutine TERMCHK is called at the end of each 


completed ICWAKE calculation. It computes and prints particular terms 
in the equations of motion, and the user can add terms as he wishes. 
TERMCHK allows the user to investigate the applicability of various 
approximations to the full equations in particular flow regions. At 
this time, TERMCHK calculates dP/dz, OR, ,/9z, dR, / er, and R,/t at 
the axis. These terms appear in the axial (W) momentum equation. 
Because of the simplicity of this subroutine, a flow chart is not 


presented here. 


SUBROUTINE N@ZERG. Subroutine NOZERO is called by DATIN to ensure that 


the turbulence quantities Ro? Ryr> Rgo> and € are nonzero (although 
small) outside of the upstream wake boundary. 

SUBROUTINE BLAYER. Subroutine BLAYER is called by ICWAKE to compute 
the profiles of the mean axial velocity and the seven turbulence quan- 
tities in the boundary layer aft of a body. Propeller effects are not 
included. See Gralidwek’, et al. (1976) for the details of the boundary 
layer algorithm. 

SUBROUTINE PROPWV. Subroutine PROPWV is called by ICWAKE to compute 
the profiles of the mean axial and swirl velocities immediately 

behind a propeller that: is attached to a body. The mean axial vel- 
ocity profile computed by subroutine BLAYER is assumed to represent 

the flow entering the propeller. See Schwartz and Bernstein (1975) 

for the details of the propeller algorithm. 

SUBROUTINE PROPU. Subroutine PR@PU is called by ICWAKE to compute the 
radial velocity profile immediately behind a propeller that is attached 
to a body. The calculation is based on the mean axial velocity profile 
compated by PR@PWV. See Grabowski, et al. (1976) for the details of 
the radial velocity calculation. 

SUBROUTINE PTURB. Subroutine PTURB is called by ICWAKE to compute the 
modification of the turbulence in the boundary layer aft of a body as 
it passes through a propeller attached to the body. The turbulence 
profiles computed by BLAYER are assumed to represent the turbulence 
coming into the propeller. See Grabowski, et al. (1976) for the details 


of the turbulence modification algorithm. 
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4. USE OF ICWAKE - General Observations 

As we pointed out in the Introduction, ICWAKE obtains a finite- 
difference solution for a rather complicated, lengthy, elliptic system 
of eleven, coupled, partial differential equations. Although this 
manual should make application of the code as easy as possible, the 
inexperienced user should be prepared for some initial difficulty. 
Several observations may prove useful. ; 

First, we note that no formal stability analysis has been 
applied to the finite-difference equations, and no theoretical 


criterion for either the maximum or optimum time-step (if such 


exists) is available. However, since the finite-difference rep- 
resentation of the equations for the mean flow and turbulence quan- 
tities is based on generally stable implicit formulations, we expect, 


and find empirically, reasonably good stability properties. In gen- 


eral, the maximum time-step that may be stably applied is at least 
two to three times the smallest mesh width in r-z space, and it is 


relatively insensitive to the free-stream Reynolds number when up- 


| 
| 


wind differencing is applied to the axial convection terms. A good 
procedure to follow in selecting a time-step is to set it to the 
smallest r-z space mesh width and, in a sequence of calculations, 
which need not be performed to convergence, gradually increase it. 
Instability caused by too large a time-step is usually catastrophic 
and readily apparent; however, we have encountered, in laminar flow 
calculations, slowly growing "wiggles" in the calculations for the 
upstream portion of the domain. These "wiggles" were smoothed in 
subsequent calculations with a smaller time-step. 
‘ The second observation is that, at sufficiently high Reynolds 
numbers when centered differences are used for the axial convection 


terms, growing or stationary "wiggles" will often appear in the down- 


stream portion of the domain in calculations. Roache (1972) discusses 
this behavior. Increasing the mesh resolution will usually rectify 
this problem, although in most cases the only practical solution is 
to switch to upwind differencing. 
A third point is that both careful consideration and experimenta- 
. tion are required in the arrangement of the computational domain and 
| finite-difference grid. We require that RMAX and ZMAX, which locate 
the NORTH and EAST boundaries, be chosen sufficiently large that 
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choosing them larger would have no substantial effect on the com- 
puted solution in the interior of the domain. At the same time we 
would like the domain to be as small as possible in order to maximize 
the possible resolution with a given, limited number of grid points. 
Thus, for a given flow or a given class of flows, some iteration is 
necessary to arrive at the optimum value of RMAX and ZMAX. 

For "production" calculations, we recommend a grid of 60 x 32 
grid points (N = 60, M= 32). For debug and trend calculations, N 
can be set to 40 and Mto 16. Optimum choice of the grid parameters | 
RCC, ZCC, NC, and MC depends very much on the nature of the flow to | 
be computed, so some physically based insight may be useful. For 
many calculations, the following values are effective: RCC = 1, 
ZCC = ZINITL + 1, MC = 8 or 14 when M = 16 or 32, respectively, and 
NC = 10 or 14 when N = 40 or 60, respectively. These values assume 
that RMAX is between about 4 and 8 and that ZMAX is between ZINITL+10 
and ZINITL+20. Flows with high swirl and large axial gradients in 


the upstream portion of the domain may require larger values of NC. 
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5. SAMPLE CALCULATION f 
The sample calculation in this section in based on the combined 
data of Swanson, et al. (1974) and Chieng, et al. (1974). These data 
were obtained in the wake behind a blunt-nosed, parallel-sided, sharp- 
sterned body with a ratio between the length and maximum diameter of 
6. The flow Reynolds number, based on the free-stream velocity and 
body radius, was 3.1 x 10°. The body was fitted with a 6-inch-diameter 
aft-mounted propeller, and the experiments were conducted under ap- 
proximately drag-free conditions. Hot-wire data were obtained at 
axial positions 4, 10, 20, 40, and 80 body radii downstream of the 
body. These include RY Rog» Ro? R and Ry,*> as well as the 
axial and circumferential mean velocity components, W and V. 
In the test case included here, we have used the experimental 
data at 4 body radii behind the body to compute the flow in the region 
behind the body to a location 20 body radii downstream. The radial ; 
boundary of tiie ~omputational domain was set at 4 body radii. Up- 
stream conditions for three flow variables were not measured in the 
experiment and were therefore estimated for this calculation as fol- 
lows: The radial mean velocity U was set to zero; the Rig turbulence 
correlation was estimated from an eddy viscosity that was obtained } 
from the known or calculable values of Roy and dW/dr; the dissipation | 
rate € was assumed to equal Kk? 72, where K = (RL + Rog + R.? 
K. + 0.53, and 2+ 0.2 (i.e., 0.2 body radii). The values of the 
constant K. and the turbulence integral scale {were suggested by 


Gran (1976). 
We wish to emphasize that the computation presented here should 


oan gr serene eveneen——yremegpeneemen 


be considered as a sample calculation, not as an accurate prediction 
of a real flow. Considerable testing and evaluation will be necessary 


before we will be able to place confidence in ICWAKE predictions. 


} 

} fe 

THE COMPUTER CODE OUTPUT INCLUDED IN THIS SECTION HAS BEEN ABRIDGED. Fr 
{ 


*This notation is defined in Grabowski (1976). 
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| Appendix | 


i The following pages contain a listing of the ICWAKE computer 
code. The routines appear in the following order: 
ICWAKE 
BLAYER 
PROPWV 
PROPU 
PTURB | 
ELINT1 j | 
. ELINT2 
QSF | 
COEFF 
NEWT | 
‘ UPCOND 
DATIN 
IUNI 
NOZERO 
. STRESS2 
PRESSR 
DIVCON 
FLXCHK 
TERMCHK 
" BLETRI 
BLKTR1 
PROD 
PRODP 
CPROD f 
CPRODP 
PPADD 
3 ® PSGF 
BSRH 
PPSGF 
PPSPF 
° COMPB 
TQLRAT 
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